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ABSTRACT 

Context. Structures seen in idealized numerical experiments on compressible magnetoconvection in an imposed strong vertical mag- 
netical field show important diff'erences with respect to those detected in observations or realistic numerical simulations of sunspot 
umbrae. 

Aims. In order to elucidate the origin of these discrepancies, here we present a series of idealized 3D compressible magnetoconvection 
experiments that diff'er from previous such experiments in several details, bringing them closer to realistic solar conditions. 
Methods. An initially vertical magnetic field Bq is imposed on a time snapshot of fully developed solar-like turbulent convection in a 
layer bounded by a stable layer from above. Upon relaxation to a statistically steady state, the structure of the flow field and magnetic 
field is examined. 

Results. Instead of the vigorous granular convection (GRC) well known to take place in nonmagnetized or weakly magnetized con- 
vection, for high values of Bq heat is transported by small-scale convection (SSC) in the form of narrow, persistent convective columns 
consisting of slender upflows accompanied by adjacent downflow patches, reminiscent of the "convectons" identified in earher semi- 
analytic models. For moderate field strengths flux separation (FXS) is observed: isolated field-free inclusions of GRC are embedded 
in a strongly magnetized plasma with SSC. Between the SSC and FXS regimes, a transitional regime (F/S) is identified where con- 
vectons dynamically evolve into multiply segmented granular inclusions and back. 

Conclusions. Our results are in somewhat closer agreement with observed umbral structures than earlier idealized models, as they do 
reproduce the strong localized, patchy downflows immediately adjacent to the narrow convective columns. Based on recent observa- 
tions of umbral dots we suggest that the conditions in sunspot umbrae correspond to the newly identified F/S transitional regime. 
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1. Introduction 

The spectacular increase in the resolution of solar observations 
experienced in the past decade has resulted in unprecedentedly 
detailed images and movies of the small scale structures seen in 
sunspots. This in turn has rekindled interest in theoretical stud- 
ies of magnetoconvection in order to understand the origin of 
these structures. This has resulted in some spectacular break- 
throughs such as the first ever numerical simulation of whole 
sunspots (Rempel et al. 2009b; 2009a). Despite such successes, 
however, there are still many questions to be answered. These 
questions include the exact nature of umbral substructure such 
as umbral dots, light bridges, dark nuclei etc., and their relation 
to structures seen in diflFerent numerical magnetoconvection ex- 
periments. 

High resolution imaging and spectropolarimetric observa- 
tions with the Swedish Solar Tower (Sobotka & Hanslmeier 
2005; Riethmiiller et al. 2008b; Sobotka & Puschmann 2009; 
Ortiz et al. 2010), Hinode (Kitai et al. 2007; Riethmiiller et al. 
2008a; Bharti et al. 2007b, 2009; Sobotka & Jurcak 2009) and 
the Dunn Solar Telescope (Bharti et al. 2007a; Rimmele 2008) 
have shown that umbral dots (UDs) have characteristic sizes of 
~ 200 km and they are mostly short-lived, the majority having 
lifetimes of a few times ten minutes or less. Larger UDs often 
split or coalesce and in high resolution images dark lanes are 
often seen inside them. In the bright parts of UDs an upflow of 
~ 1 km/s has been unambiguously detected in many instances. 
Strong localized downflow patches have been detected on the pe- 



riphery of UDs. This association with up- and downflows clearly 
demonstrates that UDs are a manifestation of magnetoconvec- 
tion in the strongly magnetized, compressible umbral plasma. 

1.1. Idealized models of magnetoconvection with an 
imposed vertical field 

In the past decades, theoretical studies of magnetoconvection 
in a strong, vertically oriented magnetic field have gone a long 
way, starting from nonlinear extensions of the classic linear sta- 
bility analysis of Chandrasekhar (1961). In the past decades, 
these studies have proceeded through the modelling of convect- 
ing magnetized Boussinesq and anelastic fluids in 2 and 3 di- 
mensions to modelling 3D compressible magnetoconvection in 
a polytropically stratified atmosphere (Weiss et al. 1996, 2002). 
One important realization in these studies was that the spatial 
separation of strongly magnetized, non-convecting plasma from 
weakly magnetized, convective flows admits overturning con- 
vection to take place also in situations where linear theory would 
predict convective stability or overstability. 

As the field strength is decreased starting from a value high 
enough to completely suppress motions, convection first sets in 
in the form of narrow, needle-like vertical convective colunms. 
The horizontal scale of these structures is small compared to the 
natural scale of nonmagnetized convection in the same layer, 
so this stage is usually referred to as "small scale convection" 
(SSC). While these convective columns can even exist as soli- 
tary "convecton" solutions (Blanchflower 1999; Blanchflower & 
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Weiss 2002; Houghton & Bushby 2010), they typically appear in 
large numbers. The "convectons" of earlier Boussinesq models 
(Blanchflower 1999; Blanchflower & Weiss 2002) are vertically 
elongated convective cells, with upflows and adjacent downflows 
embedded in stagnant fluid. In later compressible experiments 
(Weiss et al. 1996; Houghton & Bushby 2010), however, SSC 
takes the form of isolated, convective upflow plumes, their mass 
flux being compensated by a slow, uniform sinking motion in 
the surrounding magnetized fluid. The plumes are persistent and 
inmiobile, though for weaker imposed field they pulsate irregu- 
larly. 

With further reduction of the imposed magnetic field, a phe- 
nomenon known as flux separation (FXS) takes place: large iso- 
lated patches of non-magnetized, vigorously convecting fluid ap- 
pear inside the strongly magnetized component where SSC is 
still going on. The character of the vigorous convection in field- 
free patches is similar to the normal granular convection (GRC) 
in the nonmagnetic case. Flux separation was first detected in a 
simulation by Tao et al. (1998); however, it seems to be a more 
fragile mode of magnetoconvective energy transport than the 
small scale convection as in simulations without full compress- 
ibility or with lower aspect ratios it is not present. Also, mag- 
netoconvective systems show a hysteresis-like behavior where 
the presence or absence of flux separation is not uniquely deter- 
mined by the physical parameters of the system but also depend 
on its history (Weiss et al. 2002). 

For even lower field strengths, the strongly magnetized com- 
ponent becomes confined into an intermittent network of chan- 
nels separating the vigorously convecting patches. Finally, with 
the further decrease of the imposed field the topology changes, 
giving way to the well known pattern of isolated magnetic flux 
concentrations arranged in a magnetic network-like pattern be- 
tween granules (GRC). 

The above described sequence, characterized as a compe- 
tition between two modes of convection (small scale colunmal 
convection vs. vigorous granular convection), suggests an inter- 
pretation of UDs as the photospheric manifestation of SSC. At 
the same time, there are no obvious examples of flux separation 
in the solar photosphere, and the short lifetimes of observed UDs 
stand in stark contrast to the persistent character of the small- 
scale convective columns described above. 

1.2. Numerical simulations of magnetoconvection in sunspot 
umbrae 

Models like those described above are usually called "idealized" 
as they aim at a systematic study of increasingly general mag- 
netoconvective systems that are still simple enough to allow pa- 
rameter studies and to be understood in relatively simple phys- 
ical terms. In consequence, they do not aspire to faithfully re- 
produce all the physical conditions relevant for the solar pho- 
tosphere; instead, they consider convection in a poly tropically 
stratified perfect compressible fluid, with closed boundaries and 
fixed diff'usivities. 

Another approach to the problem consists in attempting to 
simulate conditions prevailing in the solar photosphere as re- 
alistically as possible, including radiative transfer and a realis- 
tic equation of state in a fully compressible equilibrium model 
(Moradi et al. 2010). An obvious drawback of this approach is 
that, going to the very limits of the available computing power, it 
is left with too few resources for systematic parametric studies. 
Yet these ambitious models have been quite successful in e.g. 
reproducing the observed penumbral fine structure. Structures 
readily identified with the observed umbral dots are also seen in 



such simulations, as first demonstrated by Schiissler & Vogler 
(2006). In contrast with the small scale convective columns of 
idealized models, these simulated UDs have finite lifetimes, and 
they show a characteristic "coff'ee bean" structure, being crossed 
by a dark lane along the long axis of their oval shape. Downflows 
are present near the end points of the dark lane. The dark lane 
itself is due to an optical depth eff'ect due to the pileup of up- 
flowing material below a magnetic cusp above the UD. Similar 
UD-like structures are also seen in the full- spot simulations of 
Rempel et al. (2009b, 2009a). 

1.3. Objectives 

Despite the impressive achievements of these state of the art 
simulations, some difl'erences between the detailed properties 
of observed and simulated umbral dots remain (cf. Bharti et al. 
2010). Even more important diflferences remain, however, be- 
tween the nature of magnetoconvective structures in idealized 
models on the one hand, and in observed or simulated sunspots 
on the other. These difl'erences include the contrast between the 
observed localized patchy downflows next to UDs vs. the lack 
thereof around small scale convective columns in experiments; 
the persistent character of convective columns vs. the limited 
lifetimes of UDs; as well as the lack of clear signs of flux sepa- 
ration in observed or simulated sunspots. 

We cannot claim a satisfactory understanding of magneto- 
convection in sunspots until the origins of these difl'erences are 
properly elucidated and a convincing link between convective 
columns and UDs is forged. For this, there is a need to bridge 
the gap between the idealized models and the realistic simula- 
tions. As a first step in this process, in this paper we construct an 
idealized magnetoconvection model where, however, we relax 
certain constraints employed in earlier idealized models such as 
Weiss et al. (2002) to bring the model closer to solar conditions. 
Specifically 

- We consider a two-layer model wherein the convectively un- 
stable fluid is bounded by a stable layer from above. 

- We use open lower boundary conditions. 

- We employ a subgrid closure scheme to account for the ef- 
fects of small-scale turbulent transport. 

- The overall stratification chosen, while politropic for sim- 
plicity, is selected to be reasonably similar to the actual solar 
case. 

- Finally, we use a novel numerical scheme, the BGK-MHD 
scheme that is particularly well suited for the capture of near- 
discontinuities in the flow like the current sheets character- 
istic of magnetoconvection. This admits a satisfactory repro- 
duction of structural details despite the relatively low resolu- 
tion used for the computation of a larger number of cases, as 
needed in parametric studies. 

In the simulations we first relax a convectively unstable layer 
of perfect gas politropically stratified under gravity. Then, a ver- 
tical uniform magnetic field is imposed on a snapshot of the 
fully developed, statistically steady turbulent convection. After 
a statistically steady state of magnetoconvection is attained, the 
properties of this system are examined. We construct a series of 
models with diff'erent strengths of the imposed magnetic field. 

Section 2 describes the setup of numerical model. Section 
3 presents and discusses the results from the simulations and 
discusses these results, while Section 4 concludes the paper. 
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2. Numerical model 

In the current study, the 3D BGK-MHD code developed by Tian 
(2010) is applied for the solution of the following resistive MHD 
equations under gravitational field: 



dpidt = -V-pt;, 

dpvldt = -V -(pvv-BB)-Vptot-^^ -^-^pg, 
dE/dt = -V • [(E + ptot)v -V Fa- BB V 

-Bxr](WxB)]-^pv g, 
dB/dt = -V • (vB -Bv- T]VB) 



(1) 
(2) 

(3) 
(4) 



where, /?tot = Pg ~^ Pm^ E = Ei -\- E^ + E^, pg is the gas pressure, 
Pm = ^B - B is the magnetic pressure, Ei is the internal energy, 
Efn = ^B - Bis the magnetic energy and E^ = \pv^ is the kinetic 
energy, g is the gravitational acceleration, Ed the diff'usive heat 
flux, E the viscous stress tensor, 77 the magnetic resistivity. All 
the other symbols have their standard meanings. The 3D BGK- 
MHD code solves the non-magnetic part of the governing equa- 
tions by BGK-NS solver modified due to the Lorentz force, the 
magnetic part by gas kinetic theory based flux-splitting method, 
using the particle distribution function constructed in the BGK- 
NS solver. The van Leer limiter is used to reconstruct the state 
around the interfaces of finite volumes. The solenoidality of the 
magnetic field (V • 5 = 0) is ensured by the field interpolation 
constraint method. The code is parallelized by the MPI protocol. 



2.1. Boundary conditions 

Side boundaries are periodic. The upper boundary is impenetra- 
ble and stress-free and the density and internal energy are fixed: 



~dz ~ ~dz ~dt ~ ~dt ~ 



(5) 



The lower boundary is transmitting, so here only the hori- 
zontal averages of density and internal energy are fixed at their 
initial values while the vertical derivatives of their fluctuations 
are set to zero: 



dv _ dp' _dp _de _ 

d~z~ ~ ~ ~dt ~ dt ~ 



(6) 



Throughout this paper, overbar and prime denote average and 
fluctuating part, i.e. for a variable a we have a = a -\- a\ Unless 
explicitly stated otherwise, overline denotes horizontal average. 
The magnetic field is constrained to become vertical 



OZ 



(7) 



at both the upper and lower boundaries. In this case no magnetic 
flux cannot be advected into or out of the computational domain 
by the flow. 

2.2. Initial setup 

Initially, the system consists two layers: an upper stable radiative 
layer (1 > z > Zi) and a lower unstable convective layer (zi > 
z > 0). The upper layer is subadiabatically stratified, the lower 
layer is slightly superadiabatic. 



nz) 



ro[i+Zi(i+z)] 
ri[i+Z2(i+z)] 



for 1 > z > Zi 
for zi > z > 



(8) 



^ ~ 1 PiiTITii^^^^^ 



for 1 > z > Zi 
for zi > z > 



(9) 



where Z\ and Z2 are the dimensionless gradients of temperature 
of the subadiabatic and superadiabatic layer, respectively. m\ and 
m2 are the polytropic indices for radiative and convective layer, 
respectively. Tq, p^ are the temperature and pressure at the top 
lid. Ti, pi are the temperature and pressure at the interface be- 
tween the initial stable and unstable layers, i.e.. 



To = T(z =1) Ti = T{zi\ 
Po = p(z = 1) Pl= p(zi). 

In hydrostatic equilibrium we have 

(mi -h l)Zi = (m2 + 1)Z2 = |^|. 



(10) 
(11) 



(12) 



The choice of Zi and Z2 determine the temperature contrast be- 
tween top and bottom, and thereby, by fitting to a standard solar 
model, the vertical extent of our computational domain. For this 
reason Z is occasionally called depth parameter in the litera- 
ture. For m2 we choose a value below the adiabatic value that 
reproduces the right value of gravity acceleration according to 
equation (12). For Z\ a lower value is adopted to make the top 
layer subadiabatic, thus convectively stable. Then mi can also 
be determined from equation (12). The specific values of these 
parameters are listed in Table 2. 

2.3. Treatment of radiation and subgrid scale motions 

Viscosity and diff'usivities due to subgrid scale turbulent eddies 
are calculated according to the Smagorinsky model. The dynam- 
ical viscosity is given by 



fi = p(c^Af(2(T : (T) 



1/2 



(13) 



where is an adjustable constant, usually chosen from the range 
0.1-0.2. The filter width A is taken to be the local resolution, 
colon stands for tensor contract and o" = diVj+djVi. The turbulent 
heat transfer coefficient Cs and the magnetic dilfusivity 77 are 
calculated as 



Cs =p/S T]= plpdn^ 



(14) 



where the Prandtl number 6 and the magnetic Prandtl number 6^ 
are assumed to be constant. 

For the current study, radiative transfer is considered in the 
diff'usion approximation. The total (radiative-Fturbulent) diff'u- 
sive heat flux is then 



Fd = -CTVT-CsVS, 



(15) 



where S = Cp(lnT - In p) is the specific entropy, Cp is the 
specific heat at constant pressure and is the adiabatic gradient. 
In the stable layer, Ct is set so that radiation carries out the input 
energy flux. In the convection zone, Ct is very close to zero. As 
Cs = p/S represents turbulent diff'usion, it is set to zero in the 
stable region. 



2.4. Dimensional analysis and model parameters 

We nondimensionalize all variables by choosing the following 
basic scales. 

The scale of the length, /sci, is the depth of the computa- 
tional domain. For the density, gas pressure and temperature, 
their respective values at the top boundary are chosen as the ba- 
sic scales, i.e., psci, Psci and T^d- The scales for the rest of the 
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Table 1. Dimensional analysis. 



Scale 


value 


unit 


Pscl 


1.03 X 10"^ 


g/cm^ 


Pscl 


9.94 X 10^ 


dyn/cm^ 


Tsci 


12992.6 


K 


hcl 


1.44 


Mm 


Uc\ 


147.22 


s 


Vscl 


9.81 


Km/s 


Bscl 


3534.8 


G 


Fscl 


9.75 X 10^1 


erg/s/cm^ 



Table 2. Summary of numerical parameters. 



Parameter 


Definition 


Value 


Grids 


N.xNyX N, 


60 X 60 X 60 


Aspect ratio 


Ix • ly • Iz 


5:5:1 


CFL number 


Ax/(v + Cs)/At 


0.3 


Depth after relaxation 


in PSHs 


3.6 


SGS Prandtl number 


S ^ jdlK 


1/3 


Deardorff number 


in Eq. (13) 


0.2 


Magnetic Prandtl number 




10 


Ratio of specific heat 


y 


5/3 


polytropic index (radiative) 


Ml 


2.728 


polytropic index (convective) 


m2 


1.485 


Depth parameter (radiative) 


Zi 


2 


Depth parameter (convective) 


Z2 


3 


Initial stable-unstable interface 




0.9 



Table 3. Imposed magnetic field strength in the numerical runs. 



Identifier 


Bo/ Be 


A 


0.35 


B 


0.70 


C 


1.34 


D 


1.41 


E 


2.12 


F 


2.82 


G 


3.10 


H 


3.53 



quantities can be deduced according to dimensional analysis. For 
instance, the scale of velocity and time are 



f^scl - (/^scl/Pscl) 



1/2 



scl 



(16) 



respectively. (Note that f;sci is the isothermal sound speed at the 
top.) In the case of magnetic field strength, an extra factor of 
In^^^ is introduced in the unit in order to give it a more physi- 
cal meaning (of thermal equipartition field strength at the top): 

Of particular importance is the scale of the input energy flux. 



scl 



Acl^scl • 



(17) 



As the solar radiation flux is fixed in physical units, its value in 
nondimensional units will depend on the value of F^d, i.e. on 
the position of the upper boundary of our box inside the solar 
convective envelope. Thus, by specifying the input flux, we can 
approximately control the radial location of our numerical com- 
putational domain within the standard solar model. 

In summary, the input energy flux, the aspect ratio and the 
depth of the domain determine the "geographic" parameters of 
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Fig. 1. Statistical properties of steady non-magnetic convection, used 
as initial condition for magnetoconvection. Shown are the vertical dis- 
tributions of (a) contributions to the total energy flux (upper left); (b) 
means (lower left) and (c) rms fluctuations (lower right) of thermody- 
namic variables; (d) rms velocity fluctuations (lower right). Averages 
are taken in time and in the horizontal plane. 
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Fig. 2. Selected histories of magnetic energy for various runs, as indi- 
cated by the labels. 



the model, i.e., size, location. The scales of density, gas pressure 
and temperature determine the state of matter. In order to per- 
form the dimensional analysis, we need a standard solar model. 
In the current study, we adopt the combined model calculated by 
Guenther et al. (1992) (Table 3B therein). 

In our thermally relaxed non-magnetic convection model, the 
input flux is ~ 0.065, implying that the top our box is at a depth 
of ~ 695 km from the solar surface (from unit continuum optical 
depth). This distance roughly corresponds to the depth of the 
Wilson depression. 

The parameters of our models are summarized in Table 2. 

For magnetoconvection it is in fact more meaningful to ex- 
press magnetic field strength in terms of the turbulent equipar- 
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tition field rather than in Alfvenic units (or thermal equiparti- 
tion field). The turbulent equipartition magnetic field Be is de- 
fined by = B^/^n. We evaluate 5^ at z = 0.91, near 
the top of the convectively unstable part but deep enough to 
be unaff'ected by the closed upper boundary condition. We find 
Be = 0.295sci = 1025 G. Table 3 presents the strength of the 
imposed vertical magnetic field in the individual runs in units of 

Be. 



3. Numerical results and discussion 

Upon starting the calculation with no magnetic field, after long 
term relaxation convection approaches a nearly steady state 
where the temporally and horizontally averaged total energy flux 
is nearly constant. The various contributions to the energy flux 
are shown in the upper left panel of Fig. 1. The strong down- 
ward directed kinetic energy flux is a well known hallmark of 
steady compressible convection, expected from convection the- 
ory (Unno & Kondo 1989; Petrovay 1990) and also consistant 
with numerical simulations (Chan & Sofia 1986). From the con- 
stancy of the total energy flux it is clear that the non-magnetic 
convection is completely relaxed. 

Next, we superpose a vertical uniform magnetic field on an 
instantaneous snapshot of the relaxed non-magnetic convection. 
Fig. 2 shows the time history of magnetic energy for some of the 
cases listed in Table 3. It is apparent that that after a short while, 
the magnetic energy contained in the system becomes statisti- 
cally steady. For the cases with very weak or very strong fields, 
the relative fluctuation of the magnetic energy is lower than in 
the intermediate cases E and F. The reasons for this will be clar- 
ified below. 

3.1. Flow morphology 

The eff'ects of a vertical magnetic field imposed on the granula- 
tion pattern are summarized in Fig. 3. 

As revealed by previous studies, e.g., Weiss et al. (1996); 
Stein & Nordlund (2000), for low field strengths magnetic fields 
are intermittent, concentrated into flux tubes confined to the nar- 
row lanes of the intergranular network (case A). For somewhat 
stronger fields the flux tubes join to form a connected network 
(case B). The suppression of motion in the area of strong mag- 
netic fields is clearly seen. 

In the opposite limit of very strong imposed magnetic field, 
convective flows are restricted to narrow columns separated by 
strongly magnetized, nearly stagnant fluid (case H). Each of 
these narrow columns consists of a single convective cell, seen 
in cross section in the figures. The size of these small scale 
convective elements is small compared to the sizes of normal 
granules, inviting comparison to the umbral dots observed in 
sunspots. This correspondence was indeed suggested by Weiss 
et al. (2002). Note, however, that the small scale convective el- 
ements in the experiments of Weiss et al. (2002) only contained 
upflows, their upwards directed mass flux being compensated by 
a slow sinking motion in the strongly magnetized bulk fluid. In 
contrast, we find strong localized downflows adjacent to the up- 
flows, within the non-magnetic convective columns. This prop- 
erty makes our small scale convective elements a close relative 
of the convectons of Blanchflower (1999) and Blanchflower & 
Weiss (2002) so for brevity, in the following we will refer to 
them by this name. 

In the intermediate cases E & F we recover the flux separa- 
tion phenomenon, first discovered by Tao et al. (1998), where 



and the nature of umbral dots 




Fig. 4. A vertical cut of our box along x = 2.5 at an arbitrary instant. The 
colour scale represents temperature; vectors represent the velocity field 
projected on the plane. Rows correspond to different runs as labeled on 
the left. 



small scale convection and vigorous granular convection coex- 
ist, the latter being present in large field free patches (in what fol- 
lows: inclusions, for brevity) embedded in a strongly magnetized 
fluid displaying SSC in the form of scattered convectons. Each 
inclusion includes several granules. The size of the inclusions is 
comparable to the size of our box: it is this efl'ect that leads to 
the large relative fluctuations of magnetic energy in cases E and 
F (Fig. 2). 

Finally, the snapshot for case G in Fig. 3 displays multiply 
segmented convective plumes of size intermediate between in- 
clusions and convectons. The lanes segmenting these elements 
are cool and coincide with downflows. This already suggests 
that our case G represents a heretofore unknown transitional 
regime of magnetoconvection between small scale convection 
(SSC) and fully developed flux separation (FXS) — a finding that 
will be corroborated below. 

Figs. 4 and 5 present vertical cuts of our box. As expected, 
stronger magnetic field is generally associated with cooler gas 
— usually even cooler than in the less magnetized downflows. 
In the weak field cases, despite our imposed boundary condi- 
tions of a vertical magnetic field, the field lines often get quite 
tangled inside the convective region as the vigorous convective 
flows freely deform them. In the strong field case, on the other 
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5 



1 




max 



min 



Fig. 3. Contour plots various physical quantities, as marked on top of each column, in a horizontal cut of our box at z = 0.78 at an arbitrary instant. 
Rows correspond to different runs as labeled on the left. 



hand, field lines remain close to vertical and the columnar shape 
of convectons is apparent. 

It is clear from these plots that the size of the various con- 
vective structures is not very small compared to the size of our 
box, leading to an imperfect cancellation of turbulent fluctua- 
tions. This results in significant temporal fluctuations of the to- 
tal magnetic energy (cf. Fig. 2). Furthermore, momentum is not 



conserved in our simulation as it can be advected in and out of 
the box through the open lower boundary. In many cases this re- 
sults in a net uniform mean flow. Both efl'ects are particularly 
strong in the flux separation cases (E and F), due to the presence 
of the inclusions with a size comparable to the domain size. In 
our plots of the velocity field the resulting mean flow has been 
subtracted for clarity. In a few runs made with larger box sizes 
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Fig. 5. Same as Fig. 4 except that vectors now represent the magnetic 
field. 



and/or closed lower boundary we confirmed that the net hori- 
zontal flow indeed vanishes in these cases, while otherwise the 
results remain unchanged. 

3.2. Time evolution 

Fig. 6 and Fig. 7 illustrate the time evolution of ver- 
tical velocity and magnetic field strength near the top 
of the convective layer. The evolution is more clearly 
seen in the corresponding animations, publicly available at 
http : //astro . elte . hu/~kris/magconvmovies 

The convectons seen in the strong field case H are fairly per- 
sistent but they slowly pulsate. The small-scale convective mo- 
tion is confined in thin vertical field- free columns. The cross sec- 
tion of these columns evolves with time periodically, as indicated 
by the small white box in first row of Fig. 6. 



It is interesting to note that as the largest convectons ap- 
proach their maximum size a downdraft often appears in their 
middle. The downdraft is obviously cooler and somewhat more 
magnetized than the rest of the convecton. This phenomenon is 
somewhat analoguous to the first phase in the development of ex- 
ploding granules as seen in quiet sun observations. In fact, this 
is but a first sign of the more spectacular temporal variation of 
the convective elements seen in case G. 

The time variation of the convective elements in case G re- 
veals a new type of behaviour apparently not seen in previ- 
ous simulations. The pulsation of some convectons reaches an 
extremely high amplitude. During their size increase, mutliple 
downdrafts appear inside them; then, these downdrafts merge 
into a network of internal downflow lanes, lending the convec- 
tive element a segmented appearance, reminiscent of a granular 
inclusion. (Again, note the analogy with exploding granules, cf. 
Hirzberger et al. 1999.) What we see here is essentially convec- 
tons developing into inclusions and back, on a time scale much 
longer than the pulsation period of normal (small) convectons. 
(In dimensionless units, the cycle time is ~ 100 vs. ~ 20.) Case 
G thus represents a hitherto unseen transitional state between 
small scale convection and flux separation. 

In the flux separation case (E and F) we have normal granu- 
lar convection going on inside the inclusions while small, fairly 
steady convectons can be observed in the magnetized fluid in 
between. There seems to be no more evolution of convectons 
to/from inclusions, though the possibility cannot be excluded 
that on time scales much longer than the length our run such 
evolution may still occur. 

Finally, for even weaker fields the magnetozed component 
becomes limited to a network of channels between the inclu- 
sions (case B), a patterns sometimes observed in plages or facu- 
lar areas. These thin magnetized channels do not support small 
scale convection any more. For the weakest fields, the magnetic 
field becomes intermittent, breaking up into flux tubes confined 
to the intergranular lanes and defining a sort of "supergranular" 
network (case A). 

3.3. Flow statistics 

For a more quantitative analysis of the properties of turbu- 
lent magnetoconvection we plot the PDFs (probability distribu- 
tion functions of the vertical components of the magnetic field 
strength (in Alfvenic units) and of the flow velocity in Fig. 8 
and Fig. 9, respectively. Only typical cases are plotted, i.e., SSC 
(case H), FXS (case E), and GRC with intermittent fields (case 
A). In Fig. 9 the non-magnetized convection case is also shown. 
The data are sampled at two horizontal layers: near the upper 
convective boundary (z ~ 0.8) and deep inside of the convective 
part (z ~ 0.33). 

The shape of PDFs for the distribution of vertical Alfvenic 
speed is similar in both layers sampled. This might be due to the 
fact that our models are not very deep and subject to vertical field 
line boundary conditions. Overall, these PDFs are rather similar 
to those obtained by Weiss et al. (2002). 

For SSC (case H, red dashed), the PDFs are relatively flat, 
with a smooth peak just below the maximum Va^ (magnetized 
component) and a sharp peak around 0.1 (convectons). In con- 
trast to Weiss et al. (2002), in the current study the range of 
field strengths extends to negative values even when the imposed 
field is strong. This may be related to the presence of localized 
downdrafts adjacent to upflow plumes in convectons: the result- 
ing strong rotational shear can play havoc with the field lines of 
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Fig. 6. Time evolution of vertical velocity (odd rows) and vertical magnetic field strength (even rows) at z = 0.8, near the top of the convective 
layer. Each pair of rows presents a diff'erent run as indicated by the labels. Cases H, G, F and E are shown. Note that in case G the time interval 
between the consecutive snapshots is significantly longer than in the other cases. For case H a white box encloses a pulsating convective tube. 
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Fig. 7. Same as Fig. 6, for cases D, B and A 



the weak magnetic field inside the convective column, occasion- 
ally turning around their vertical orientation. 

For the intermittent field case with GRC (case B, black solid) 
the form of the PDF reflects turbulent flux expulsion, the sym- 
metric peak around the origin corresponding to the weakly mag- 
netized, vigorously convecting fluid while the smooth secondary 
peak at higher field strengths is due to the magnetic flux concen- 
trations. This form of PDF is familiar both from previous simu- 
lations and from quiet sun observations (de Wijn et al. 2009). 

For FXS, (case E, magenta dotted) the shape of the PDF is 
similar to case B but the strong field peak is higher and shifted 
towards higher field strengths, as expected. 

Unlike the PDFs of the magnetic field strength, the PDFs 
of the flow velocity sampled in difl'erent layers show obvious 
discrepancies, as obvious from Fig. 9. 



For GRC with intermittent magnetic field (case B, thick 
black solid), the PDFs of vertical velocity are almost identical 
to those of non-magnetized convction (thin blue solid). Their 
change of shape between the top and the middle of the convect- 
ing layer is related to the well known change of morphology 
from the observed granular pattern near the surface to isolated 
fast downdrafts embedded in a slowly rising, warm bulk fluid in 
the bulk of the convective zone. The lack of a clear secondary 
peak at negative velocities in the top panel may be due to the 
low filling factor of the downdrafts and possibly also to a lack of 
a characteristic scale among them. 

As the imposed magnetic field is increased, the main peak of 
the PDFs shifts below zero. This indicates a general slow sink- 
ing of the magnetized fluid component. An extended shoulder 
on this peak on the negative side corresponds to the faster down- 
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Fig. 8. PDFs for the distribution of vertical Alfvenic speed, Va^ = 
Bz/ yj^7[p, at z = 0.8) and z ~ 0.33 for various imposed vertical mag- 
netic field strengths as indicated inside the labels, corresponding to SSC 
(red dashed), FXS (magenta dotted) and GRC (black solid). 

drafts localized next to the upflows. Both this shoulder and the 
secondary peak on the positive side representing the upflows are 
more distinct near the top of the convecting layer. 

The anisotropy of turbulence is shown in Fig. 10. As ex- 
pected, the ratio of vertical component to total field strength is 
approaching 1 as we increase the imposed magnetic fields (right 
hand panel). Even for the very weak field case, this ratio never 
goes below 0.9. This is simply because the field lines can only 
be significantly deflected from the imposed vertical direction in 
places where the magnetic field is inherently weak. 

On the left hand panel the curves can be naturally grouped 
into three classes, according to their shape and position. These 
classes correspond to SSC (case H), FXS (E and F) and GRC 
with intermittent field (cases A-D). This lends more objective 
support for considering SSC, FXS and GRC three physically 
diflTerent modes of magnetoconvaction. Furthermore, the curve 
representing case G (red dashes) is quite obviously a "hybrid", 
its shape resembling case H but its overall position placing it in 
the FXS group. This, in turn, lends further support to our iden- 
tification of case G as a new transitional mode between the SSC 
and FXS modes of magnetoconvection. 

3.4. Comparison with observations 

In the weakly magnetic cases the morphology, time dependence 
and statistics of our magnetoconvective flow is in good accor- 



Fig.9. PDFs for the distribution of vertical velocity, i.e.,v^, at z = 0.8) 
and z ~ 0.33 for various imposed vertical magnetic field strengths as 
indicated inside the labels, corresponding to SSC (red dashed), FXS 
(magenta dotted) and GRC (thick black solid). The blue thin solid thin 
is the PDF of the non-magnetized case. 



dance with observations. This is so despite the fact that, when 
scaled back to physical units along the lines of Section 2.4, our 
convective structures are significantly larger and longer lived 
than those observed in the solar photosphere. We attribute this 
diflference mainly to the lack of a freely radiating photosphere 
in our experiments, which is known to have a key role in de- 
termining granular scales (Rast 2003). The overall qualitative 
and quantitative similarity of the magnetoconvective structures 
in our weak field models to observations encourages us to look 
for similar observational parallels with the structures seen in our 
strong field experiments. 

The analysis of high-resolution Doppergrams and continuum 
images of sunspot umbrae by Bharti et al. (2007a) show that the 
umbral dots are surrounded by downwflows. Localized down- 
flows near UDs have also been observed by Ortiz et al. (2010). 
These observed umbral structures are quite similar in appearance 
to the convectons surrounded by a patchy ring of downflows, 
seen in our experiments in the SSC domain. 

Furthermore, in larger umbral dots dark lanes are often ob- 
served. These UDs rarely have the regular "coflTee bean" struc- 
ture suggested by the simulations of Rempel et al. (2009b); 
instead, dark lanes are often asymmetrically positioned inside 
them, and multiple dark lanes are common, lending a segmented 
appearance to these UDs. Umbral dots with up to six dark 
lanes have been observed by Bharti et al. (2007b). Observations 
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Fig. 10. Vertical anisotropy of the velocity field (left hand panel) and magnetic field (in Alfvenic units, right hand panel), for cases A (thin black 
solid), B (green dotted), D (green dashed), E (yellow dash-dot), F (orange dash-3dot), G (red long dashes) and H (red solid). 



also indicate that the UDs evolve rapidly compared to sunspots. 
Bharti et al. (2007b) showed how a UD with multiple dark lanes 
changed its appearance over the course of about an hour. (Note, 
however, that some of the apparent temporal variations in UDs 
may be due to seeing eflTects, cf. Hamedivafa 201 1.) 

These latter properties of UDs are more reminiscent of the 
transitional convective structures seen in our case G. Indeed, 
case G represents a vertical magnetic field strenght of about 3 kG 
in physical units, just in the range expected for sunspot umbrae. 
We therefore suggest that conditions in sunspot umbras corre- 
spond to the FXS-SSC (F/S for short) transitional regime identi- 
fied in this paper. 

4. Conclusion 

We have constructed a series of numerical experiments of mag- 
netoconvection in a layer of compressible plasma with an open 
lower boundary and bounded by a stable layer from above, us- 
ing a BGK-MHD scheme. Our model setup diflTered from that 
of previous idealized magnetoconvection models (Weiss et al. 
2002) in several important aspects, being closer to realistic solar 
conditions. Despite this, the general behaviour of the solutions 
as a function of the imposed field strength does basically agree 
with those models. In particular, flux separation (FXS) still pre- 
vails in a wide parameter range and the small scale convection 
(SSC) dominates the energy transport in the strong field case. 

An important diflTerence compared to the results of Weiss 
et al. (2002) is that in the SSC regime strong, narrow, patchy 
downflows appear in a ring around the narrow upflows in our ex- 
periments. Thus, our "convective columns" may be more closely 
related to the solitary "convectons" of Blanchflower & Weiss 
(2002) than to the localized plumes in previous magnetoconvec- 
tion experiments. This property of our experiments is in agree- 



ment with recent high resolution observations of sunspot umbrae 
(Bharti et al. 2007a; Bharti et al. 2009; Ortiz et al. 2010). 

Furthermore, we have identified a hitherto unknown transi- 
tional regime between FXS and SSC. This F/S regime is very 
dynamic, with a continuous evolution of convectons into gran- 
ular inclusions and back. Intermediate states of the convective 
elements are often reminiscent of larger, multisegmented um- 
bral dots, also known to be subject to dynamic changes. This 
suggests that conditions in sunspot umbrae correspond to the 
F/S transitional regime identified in this paper and that UDs es- 
sentially represent convectons and/or convecton/inclusion tran- 
sitional phases as seen in our run G. 

It should be noted that the realistic numerical simulations of 
Schiissler & Vogler (2006) and Rempel et al. (2009b) suggest 
a somewhat diflTerent interpretation of the substructure and time 
dependence seen in UDs. The dark lanes in these simulations 
are due to an optical depth eflfect and they do not coincide with 
downflows; downflows, however occur near the end points of 
the lanes. The finite lifetime of UDs, in turn is due to the "safety 
valve" mechanism operating at the magnetic cusp above the con- 
vecting columns, in a manner somewhat akin to the suggestion 
of Choudhuri (1986). High resolution observational studies of 
the correlation and relative position of downflows and dark lanes 
may help to clarify this issue in the future. 

It should be added that while observations suggest that the 
properties of central umbral dots are diflferent from those of pe- 
ripheric umbral dots (Sobotka & Jurcak 2009), it is hard to draw 
a sharp line between the two groups and the finite tilt of the mag- 
netic field may have an important role in determining the proper- 
ties of central UDs as well. Non- vertical fields, e.g. due to a fan- 
ning out or twisting of the flux rope, may also be instrumental in 
inhibiting the appearance of flux separation inside sunspot um- 
brae — though this may be also explained simply by the strong 
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fields in sunspot umbrae or by the large depth (and consequently 
low aspect ratio) of the umbra. 

These considerations suggest that one possible way to bring 
idealized models closer to realistic numerical simulations could 
be the extension of these results to the case when the imposed 
magnetic field is tilted. 

Acknowledgements. The computations were tested and performed on Eotvos 
University's 416-core, 3.7 Tflop HPC cluster Atlasz and on the supercomput- 
ing network of the NIIF Hungarian National Supercomputing Center (project 
ID: 1117 fragment). Support by the Hungarian Science Research Fund (OTKA 
grant no. K67746), by the European Commission's 6th Framework Programme 
(SOLAIRE Network, MTRN-CT-2006-035484) as well as by the European 
Union with the co-financing of the European Social Fund (grant no. TAMOP- 
4.2. 1/B- 09/1/KMR-2010-0003) is gratefully acknowledged. 



References 

Bharti, L., Beeck, B., & Schiissler, M. 2010, A&A, 510, A12+ 

Bharti, L., Jain, R., & Jaaff'rey, S. N. A. 2007a, ApJ, 665, L79 

Bharti, L., Joshi, C, & Jaaff'rey, S. N. A. 2007b, ApJ, 669, L57 

Bharti, L., Joshi, C, Jaaff'rey, S. N. A., & Jain, R. 2009, MNRAS, 393, 65 

Blanchflower, S. 1999, Physics Letters A, 261, 74 

Blanchflower, S. & Weiss, N. 2002, Physics Letters A, 294, 297 

Chan, K. L. & Sofia, S. 1986, ApJ, 307, 222 

Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford: 

Clarendon Press) 
Choudhuri, A. R. 1986, ApJ, 302, 809 

de Wijn, A. G., Stenflo, J. O., Solanki, S. K., & Tsuneta, S. 2009, Space Sci. Rev., 
144, 275 

Guenther, D. B., Demarque, P., Kim, Y., & Pinsonneault, M. H. 1992, ApJ, 387, 
372 

Hamedivafa, H. 2011, Sol. Phys., 45 

Hirzberger, J., Bonet, J. A., Vazquez, M., & Hanslmeier, A. 1999, ApJ, 527, 405 

Houghton, S. M. & Bushby, R J. 2010, ArXiv e-prints 

Kitai, R., Watanabe, H., Nakamura, T., et al. 2007, PASJ, 59, 585 

Moradi, H., Baldner, C, Birch, A. C, et al. 2010, Sol. Phys., 267, 1 

Ortiz, A., Bellot Rubio, L. R., & Rouppe van der Voort, L. 2010, ApJ, 713, 1282 

Petrovay K. 1990, ApJ, 362, 722 

Rast, M. R 2003, ApJ, 597, 1200 

Rempel, M., Schussler, M., Cameron, R. H., & Knolker, M. 2009a, Science, 325, 
171 

Rempel, M., Schussler, M., & Knolker, M. 2009b, ApJ, 691, 640 
Riethmuller, T. L., Solanki, S. K., & Lagg, A. 2008a, ApJ, 678, L157 
Riethmuller, T. L., Solanki, S. K., Zakharov, V., & Gandorfer, A. 2008b, A&A, 

492, 233 
Rimmele, T. 2008, ApJ, 672, 684 
Schussler, M. & Vogler, A. 2006, ApJ, 641, L73 
Sobotka, M. & Hanslmeier, A. 2005, A&A, 442, 323 
Sobotka, M. & Jurcak, J. 2009, ApJ, 694, 1080 
Sobotka, M. & Puschmann, K. G. 2009, A&A, 504, 575 
Stein, R. F & Nordlund, A. 2000, Sol. Phys., 192, 91 
Tao, L., Proctor, M. R. E., & Weiss, N. O. 1998, MNRAS, 300, 907 
Tian, C. 2010, ArXiv e-prints=1010.3594 
Unno, W. & Kondo, M. 1989, PASJ, 41, 197 

Weiss, N. O., Brownjohn, D. P, Matthews, P C, & Proctor, M. R. E. 1996, 

MNRAS, 283, 1153 
Weiss, N. O., Proctor, M. R. E., & Brownjohn, D. R 2002, MNRAS, 337, 293 



12 



